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ABSTRACT 

The energy diffusion coefficients D„(E) = {AE(E)" / At) (n = 1,2) for a system of 
equal mass particles moving self-consistently in an #-body realisation of a King model 
are computed from the probability per unit time, V(E, AE)dEdAE , that a star with 
initial energy E will undergo an energy change AE. In turn, V is computed from the 
number of times during the simulation that a particle in a 'state' of given energy un- 
dergoes a transition to another state. These particle 'states' are defined directly from 
the time evolution of E by identifying them with the event occuring between two local 
maxima in the E(t) curve. If one assumes next that energy changes are uncorrelated be- 
tween different states, one can use diffusion theory to compute D„(E). The simulations 
employ N = 512, 2048, • • • , 32768 particles and are performed using an implementation 
of Aarseth's direct integrator #-bodyl on a massively parallel computer. The more 
than seven million transitions measured in the largest N simulation provide excellent 
statistics. The numerically determined _D(£')'s are compared against their theoretical 
counterparts which are computed from phase-space averaged rates of energy change due 
to independent binary encounters. The overall agreement between them is impressive 
over most of the energy range, notwithstanding the very different type of approximations 
involved, giving considerable support to the valid usage of these theoretical expressions 
to simulate dynamical evolution in Fokker-Planck type calculations. Even so, diffusion, 
as judged from these measurements of the diffusion constants, is stronger than expected 
from theory, both in core and outer halo, by a factor up to two, rather independent of 
particle number. The experimental _D's obey very well the expected scaling cx #/lnA 
with particle number N . 

Key words: stellar dynamics - globular clusters: general - open clusters and associ- 
ations: general - numerical methods - celestial dynamics 



in 



> 
o 



1^ 

O 




in 

On 

Oh 
I 

O 



1 INTRODUCTION 



The energy per unit mass E = l2 — 4> of a, star moving 
with velocity v through a stellar system is not conserved 
when the potential 4> changes in time. Temporal changes in 
4> are due to the whole spectrum of modes in the system 
which range from collective modes (which also occur in the 
limit N oo, e.g. Weinberg 1993), over multiple encounters 
all the way down to binary encounters. This redistribution 
of energy over the different particles may lead to structural 
changes of the system. 

Following Spitzer (1987), the dominant contribution to 
the changes in particle energy for a system in dynamical 
equilibrium comes from the cumulative effects of many dis- 
tant binary encounters, each changing i? by a small amount 
and causing E{t) to perform a random walk in energy space 
(see Spitzer 1987, p. 29 for a discussion). In this case it is 
appropriate to describe the rate of change of N{E, t)dE, the 
number of particles with energy between E and E + dE, by 
a diffusion equation. 



A theoretical expression for the diffusion coefficients in 
a system of particles with forces and characterised by 

an arbitrary distribution function is given by Rosenbluth et 
al. (1957). The derivation assumes two-body encounters are 
independent and sums the contribution from encounters of 
given impact parameter b and impact velocity v. Unfortu- 
nately, the resultant expressions diverge logarithmically in 
the gravitational case due to the contribution of distant en- 
counters because, unlike in the plasma case, such distant en- 
counters are not screened. These diffusion coefficients- with 
an imposed cut-off, the Coulomb logarithm- are widely used 
in Fokker-Planck type calculations (see e.g. Chernoff and 
Weinberg (1990) or Spitzer (1987) for references). 

In a real system, a particle is undergoing several en- 
counters with a range of impact parameters simultaneously. 
Consequently, the approximation where one just adds the 
contribution from all impact parameters for pure two-body 
motion is questionable. Indeed, for the dominant encoun- 
ters with intermediate impact parameter, the description of 
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the motion as being a pure two-body encounter is likely to 
be an oversimplification. In addition, the description likely 
overestimates the importance of more distant encounters, 
both in homogeneous and in concentrated systems. In a ho- 
mogeneous system for example, it would seem that the many 
encounters with impact parameter h ^ I (for some suitably 
chosen I) perturb the orbit of a particle under consideration 
to such an extent that encounters with h ^ I are effec- 
tively terminated, thereby quenching the contribution from 
encounters with such larger impact parameters. This argu- 
ment is even stronger in a concentrated system: consider 
encounters with large h between a core and a halo particle. 
The dynamical time scale of the core particle tdyn ~ 
with p the density at the position of the core particle, may 
be much smaller than the typical encounter time tenc ~ i/v, 
with V the relative velocity. Again, this leads to a quenching 
of the contribution of encounters with large impact param- 
eter. In practice, ignorance of the appropriate value of the 
impact parameter of encounters that still contribute to the 
diffusion process is hidden in the value of the Coulomb log- 
arithm used. 

Similar considerations led Chandrasekhar (1941) to con- 
sider a different formulation of the problem based on the 
notion of particle 'states' . A particle is subject to the fluc- 
tuating gravitational field due to all the other particles caus- 
ing E{t) to perform a random walk. The force F{t) acting 
on a particle is correlated over short times but uncorrelated 
over longer times, due to the chaotic nature of the system. 
A particle is considered to be in a particular state as long 
as F is strongly correlated in time. It may then undergo a 
transition to another state and F will not be correlated be- 
tween these different states. The rate of diffusion in energy 
space of the particle can now be computed from knowledge 
of the average lifetimes of states T{F) and the probability 
distribution W{F) of a given force acting during a state. 
Maybe somewhat surprisingly, the resultant form of the ex- 
pression for the diffusion coefficient is identical to that of the 
other formulation based on independent binary encounters, 
and their numerical values differ only by a factor Si 1.11 
(Chandrasekhar 1941). 

The aim of this paper is to evaluate the diffusion co- 
efficients D„ (n = 1,2) as a function of E from a direct 
TV-body simulation by studying the properties of the ran- 
dom walk in energy space for particles of given energy in a 
similar spirit as Chandrasekhar's theory of states. We de- 
fine the states of the particle from properties of the E{t) 
curve and gather statistics on the lifetime and transition 
probability of states to compute V{E, AE), the probabil- 
ity per unit time that a particle of energy E undergoes a 
transition to another state of energy E + AE. Comparison 
can then be made between the numerically determined D's 
and their theoretical counterparts as used in Fokker-Planck 
calculations. The advantage of the numerical method over 
the analytical derivation is that the full non-linear dynam- 
ics is treated consistently, properly taken into account the 
quenching of distant encounters, the effect of interactions 
during simultaneous encounters and the possible contribu- 
tion from collective modes. 

The diffusion coefficient 1)2 (-E) can be used to define 
the relaxation time Te{E) of the system (Eq. (20) below). 
Previous measurements of Te from TV-body simulations have 
used a variety of other methods (see Huang et al. 1993 for 



more details), e.g. the rate of energy exchange between dif- 
ferent mass components, the mean-squared energy change 
of equal mass particles and the measurement of deflections 
of test stars moving through N field stars. 

The rest of this paper is organised as follows: after defin- 
ing the diffusion coefficients and giving their standard theo- 
retical expressions in section 2, we describe the experimen- 
tal set-up (3) and the experimental definition of the D's (4). 
These are compared in section 5. The relative merit of the 
proposed method for measuring diffusion coefficients and the 
reason for the (small) discrepancies between theoretical and 
numerical measurements is discussed. Finally, the paper is 
summarised. 



2 THEORETICAL DESCRIPTION 

2.1 Theoretical energy diffusion coefficients 

Spitzer (1987 and references therein) gives the standard ex- 
pressions for the diffusion coefficients based on assuming in- 
dependent binary encounters. Let N{E,t)dE be the number 
of stars in a system with energy in the range E, E + dE at 
time t. The encounter term in the Fokker-Planck equation 
is written as (Spitzer 1987, Eq. (2-71)) 



dN(E,t) 
dt 



d 



{N{E,t)D,{E)) 



dE 

1 ^ 

"2 dE^ 



{N(E,i)D2(E). 



(1) 



The diffusion coefficients D can then be computed from 
phase-space weighting the energy changes (AE" / At) (n = 
1, 2) undergone by a star per unit time. The phase-space 
weighted average {A)v of a quantity A is (Spitzer 1987, 
Eq. (2-80)): {A)v = Jo'""" A vr^dr/p, with p, the phase- 
space volume accessible per unit interval of E, given by 
Eq. (4) below. The expressions for (AE" / At) are given by 
Eqs. (2-50, 2-51) in Spitzer (1987). (We write {AE / At) in- 
stead of Spitzer's {AE) to stress these are energy changes 
per unit time.) If the background stars have masses m/ then 
the diffusion coefficients for a star of mass m are given by: 



Di(E) 

and 

D2{E) 

Here 
p{E) 



167r^m|^lnA( / f{Ef) dEf 



m fp{E) 



f{Ef ) p{Ef ) dEf 



(2) 



327r^m^ In A ' 



p{E) 



f{Ef) q{Ef) dEf 



(2(£ -*))'/" dr 



ax(-B) 



(2(£ -*))'/" dr, 



(3) 

(4) 
(5) 
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and rmax(-E) is the maximum radial distance a star of energy 
E can wander and In A is the 'Coulomb logarithm' . One has 
In A Si Inpmax/pmin whcrc pmax and pmin dcuotc maximum 
and minimum impact parameter for binary encounters, re- 
spectively. The distribution function / occurring in these 
expressions is normalised such that its integral over phase- 
space equals N. The collision term of the Fokker-Planck 
equation (2-86) in Spitzer (1987) is obtained upon substi- 
tution of Eqs. (2) and (3) into (1). The diffusion coefficients 
for a system of equal masses are given by: 



ME) 



and 



D2{E) 



IGTr^MMn A / 1 



N 



N 



f(Ef) dEf 



Np{E) 



f{Ef)p{Ef) dEf 



(6) 



327r^M^ln A / 



N 
Np{E) 



\Np{E) 

CO 

fiEf) dEf 



f(Ef) q(Ef) dEf 



(7) 



with M = Nm. 



2.2 Coulomb logarithm 

The Coulomb logarithm In A occurs because of the diver- 
gence of the diffusion coefficients due to the cumulative ef- 
fects of many distant encounters. In comparing the numeri- 
cal D's, defined in the section 4.3 below, with the theoretical 
expressions given previously, we will take pmax = Tc (Spitzer 
1987, p. 28), with Tc the core radius of the model, although 
other choices could be made as well (e.g. Spitzer (1987, p. 30) 
and Farouki and Salpeter (1994) argue for pmax equal to 
the half mass radius). In a system without smoothed grav- 
itational forces one usually takes po, the impact parameter 
causing a 90° deflection (Eq. (16)), for the minimum impact 
parameter. This gives: 



lnA = ln^:i^Ml^.ln(A.n 

2M V i 7. 



(8) 



with f (0)^ characterising the central velocity dispersion. 

In the numerical simulations presented here, gravita- 
tional forces are smoothed with e = 0.5rf (see Eq. (13) be- 
low) with d, the average central interparticle distance, given 
by Eq. (14) below. Substituting e for the smallest impact pa- 
rameter pmin gives the Coulomb logarithm appropriate for 
the numerical simulations: 



In A = In — = ln(A2Ar^/^). 



(9) 



Numerical values for Ai and A2 applicable to the King mod- 
els studied in section 3.1 are: 

Ai = 0.028, A2 = 0.019 for PFo = 9 and (10) 
Ai = 0.196, A2 = 0.197 for Wo = 3. (11) 

3 EXPERIMENTAL SET-UP 
3.1 A^-body model 

All results presented here come from simulations of A^-body 
realisations of King models (e.g. Binney and Tremaine 1987, 



henceforth BT87, p. 232). King models form a one parame- 
ter family of 'lowered isothermal' models whose distribution 
function / has a sharp cut-off at the 'tidal' energy Eo < 0. 
They can be characterised by the ratio Wo = '^{0)/a^ of 
the central potential ^'(0) over a parameter a characteris- 
ing the velocity dispersion. Given the run of density, po- 
tential and velocity dispersion for a given Wo , a particu- 
lar A^-body realisation of this King model is made using a 
random number generator. For small particle numbers, the 
resultant system may be slightly out of equilibrium due to 
small A^ statistics. In addition, the equations of motion are 
integrated by softening the gravitational force (section 3.3) 
which causes the initial state of all A^-body realisations to 
be slightly out of equilibrium. The dynamical (or crossing) 
time used in the following is defined in the usual way as 
td = tcr = M^/^/{2\ET\f'^, with M and Et the total mass 
and total energy of the cluster, respectively. (In the 'stan- 
dard' Af-body units, where M = G = -AE = 1, td = 2^2.) 
Here and in the following we take the gravitational constant 
G = 1. 



3.2 Af-body code 

The Newtonian equations of motion were integrated using 
Aarseth's A^-bodyl code (Aarseth 1985) in 15 digit precision 
using force smoothing and not including regularization. A^- 
bodyl is a high-order scheme which uses Newton divided 
differences to compute e.g. the force F,{t) on particle i as a 
function of time t as a Taylor expansion in t including terms 
up to the fourth order. Such high-order expansions are also 
used to update the positions and velocities of particles using 
individual time steps 8t, , which are computed from (Aarseth 
1985, Eq. (9)): 



|ir||ir(2)| + |j^(i)|2 



(12) 



\F(i)\\F(^)\ + \F(^)\^' 

where F^'^ = d'F{t)/dt' and rj = 0.03 is an accuracy param- 
eter. 

This code was implemented on a massively parallel 
computer, the 8192 processor Connection Machine 2 at the 
Scuola Normale Superiore di Pisa. Parallelism was exploited 
in computing the force F, on particle i due to all other parti- 
cles j in parallel using the FORTRAN 90 intrinsic function 
SUM. A catalog scheme was used to group particles with 
time steps equal to within a factor of two in bins. Such 
binning allows parallelism in updating the Newton divided 
differences and in addition significantly improves interpro- 
cessor communication efficiency. The resulting code runs at 
near 380 Mflops, about 25% of the theoretical peak perfor- 
mance, with typically Si 97% of CPU time spent in pure 
force evaluation (see Theuns and Rathsack 1993 for more 
details). 

3.3 Force calculation 

The force F,^ on particle i due to particle j is computed 
using direct summation and is softened according to: 



F,^ 



(13) 



since otherwise large integration errors occur in the absence 



of regularization, due to the singularity at 



The 
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size of the applied smoothing can be compared with several 
scales in the cluster: the average interparticle distance at 
the centre d, the semi-major axis a of a binary just on the 
division line between being soft and being hard, (e.g. BT87 
p. 534, Phard being the period of this binary), and the impact 
parameter that causes a 90° deflection, po (Spitzer f987, 
Eq. (2-5)): 

d ={mlp{Q)fl\ (f4) 
a = m/2v^, (f5) 
Po = 2m/vl^, (f6) 

where p{0) is the central density, Vm the local velocity dis- 
persion and Vco the velocity at infinity for the two-body 
encounter leading to a 90° deflection. The dependence on N 
can be brought out by defining a linear dimension R of the 
model from p(0) = M/{^R^) and by writing vl^ x 2M/R 
and in addition substituting Vco by Vm- This translates 
Eqs. (f4-f6) into d = {-i-K Izf'^ R/ N^l'-" and a = po/4 = 
R/4:N. In the simulations presented, e = d/2 oc N^^^po, 
from which it is clear that strong encounters are completely 
suppressed for large N. 



4 EXPERIMENTAL DESCRIPTION 
4.1 Basic assumption 

In the following we consider a given particle to remain in the 
same state in the time interval between two local maxima 
of the E{t) curve of that particle. We identify the energy E 
at the first of the two maxima as the state's 'energy' and 
the duration 8t of the interval as the state's lifetime. At the 
second maximum the particle makes a transition to a new 
state with in general different energy and lifetime. In the 
following we show how to compute the D's from knowledge 
of the state's energies, lifetimes and transition energies. 

This heuristic definition of a state implicitly assumes 
that energy fluctuations cease to be correlated between suc- 
cessive energy maxima, due to the chaotic nature of the 
TV-body orbits. The typical lifetimes of the states defined 
in this way turn out to be of the order of the ratio of the 
interparticle distance over the local velocity dispersion, in 
line with Chandrasekhar's (1941) estimate. 

The prime advantage of this definition is that the sam- 
pling rate used to determine the D's is adapted automati- 
cally to the particle's rate of energy change: a particle un- 
dergoing frequent energy changes (a core particle) will be 
sampled at a higher rate than a particle undergoing fewer 
energy changes (a halo particle), thereby greatly improving 
the statistics on transitions between states. In addition, the 
expected systematic energy change, due to the diffusion pro- 
cess itself, is negligible with respect to the particle's energy 
over the lifetime of the state, so that the diffusion process 
can actually be studied as a function of the particle's energy. 
Finally, this method allows one to measure the diffusion rate 
for particles in a system over a time span very small com- 
pared to the actual relaxation time of that system because 
it is based on measuring statistics of transitions and not ac- 
tual changes in structure. In practice it suffices to gather 
transition statistics over several dynamical times, to allow 
state lifetimes to be properly sampled. This enables one to 
employ a much larger number of particles in the simulations 



and numbers realistic for a globular cluster are well within 
reach of present day supercomputers. 

However, it is unclear to what extent the assumption 
of successive states being uncorrelated is satisfied. McMil- 
lan et al. (1988) also studied the properties of the random 
walk to derive a diffusion coefficient. They used the Fourier 
transform of the energy time series E{ti), where single par- 
ticle energies E were sampled at times = iSt, i = 1, ■ ■ ■ M . 
An advantage of their method is that they can tailor the 
sampling interval via its Fourier transform to establish that 
the diffusion limit is reached, i.e., that the long time be- 
haviour AE^ oc t characteristic for diffusion is sampled. A 
disadvantage of their method is that they only measure the 
second diffusion coefficient D2 and in addition have some un- 
certainties about which energy the measured D2 should be 
associated with, since the energy of the particle may change 
significantly between to and tM- 

In the following we will consider the identification of 
states as an Ansatz to be born out by further analysis. 

The method described here gives wrong results for par- 
ticles in a bound binary system*. The argument goes as 
follows: a member of a binary will not have constant energy 
(unless the eccentricity is zero). Consequently, the method 
described so far will wrongfully decide that such a particle is 
undergoing frequent transitions (one per period) but clearly 
energy changes are correlated from one state to the next. 
Fortunately, dynamically formed binaries should not be a 
problem in this investigation because they are soft due to 
the applied softening and hence will not strongly influence 
measurements of the D's. 

4.2 The transition probability V{E, AE) 

The energy per unit mass E{t) and its first three derivatives 
as a function of time t are computed during the simulations 
for all particles using the same high-order method as used 
to update positions and velocities. Local energy maxima are 
detected from dE/dt = 0, d^ E / dt^ < and are used to de- 
fine beginning to and end ti of particle states. Elapsed time 
8t = ti — to, initial particle energy E{to) and the energy 
change AE{%) = 100 * {E{ti) - E{to))/E{to) are recorded 
and stored for each transition. This is done by counting the 
number of transitions cr{E, AE, 8t) with energy change be- 
tween AE and AE + dAE, from a given state characterised 
by an initial energy between E and E + dE and a lifetime 
between 8t and 8t + d8t. Bin boundaries for storing this num- 
ber of transitions a are chosen as follows: initial energy E is 
binned linearly from i?min = -1.5*(0) to i?max = 0.01*(0), 
relative energy change AE / E is binned logarithmically from 
0.01% to 50%, and state lifetimes are binned logarithmically 
from 0.1 Phard to 3td. We use 128 bins in E,2x 128 -|-1 bins 
for AE I E (128 for positive energy changes, 128 for negative 
energy changes and 1 bin for \AEI E\ < 0.01%) and 128 bins 
for 8t. Given a, V can be computed from 

N(E)V(E, AE)dEdAE = <y{E, AE, 8t), (17) 

St 

where dE and dAE are bin widths and T is the total simu- 
lation time. 

* I would like to thank S. Aarseth for pointing this out to me. 
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Table 1 . Summary of runs 

N/512 Wo T/ti " ^^EtIEt'' 



1 

4 
4 
16 
64 



9 4 X (10+4) 

9 27+3 = 

3 27+3 = 

9 9 

9 9 



8.2E-5 127 855 

1.8E-3 504 868 

1.2E-3 297 687 

1.4E-3 1 075 069 

2.0E-3 7 154 739 



^ total simulation time in dynamical times 

^ Relative total energy change over simulation time 

^ Total number of transitions during run 

^ four runs averaged, each system was run over 4 before 
collecting statistics 

^ integrated over 3 before collecting statistics 



4.3 Experimental energy diffusion coefficients 

Several quantities of interest can be computed as moments 
of a, or alternatively of V, through Eq.(17). The number 
N{E)dE of particles at the energy E can be computed from: 



N(E)dE = 7^ ^ (^(E, AE, 8t) 8t, 



(18) 



where dE is the energy bin width. In turn, the diffusion 
coefficients can be computed from 



Dn{E)N{E)dE = ^ ^ a{E,AE,8t) (AE)" 



(19) 



Note that these equations are only accurate if the major 
contribution to the sums on the rhs is due to states with 
lifetimes 8t < 5tmax = 3td << T, since only those are sam- 
pled properly. 



5 RESULTS 

Table 1 presents a summary of the runs performed. The CPU 
time required to perform these runs is dominated by the 32k 
model which needed approximately 10 days of CPU time to 
complete. The time evolution of the Lagrangian radii was 
investigated to look for obvious signs of evolution in these 
models. Any systematic changes in these Lagrangian radii 
are both small with respect to the erratic changes due to 
low particle number and with respect to the difference in 
Lagrangian radii between the Wo = 3 and Wo = 9 models. 

Figure 1 compares the number of particles N{E)/N ob- 
tained from Eq. (18) against the corresponding theoretical 
curve for the initial King model. The agreement between 
them is excellent which shows that time averages of the sys- 
tem obtained from the simulations are a good measure for 
properties of the initial state of the system. We suggest that 
any small differences are due to initial transients, caused by 
the fact that the initial models are not completely in equi- 
librium. 

A comparison of the theoretical diffusion rates 
N(E)Dn(E) (n = 1,2) with the D's from Eqs. (6) and 
(7) respectively, using the value for the Coulomb logarithm 
computed from Eq. (11) against the experimental values as 
defined in Eq. (19) is made in Fig. 2 for the N=32k, Wo = 9 
case (but see also section 6.3 on mean field relaxation). The 



dimensionless D's in this figure are scaled so as to make 
them independent of N. The overall agreement is impres- 
sive in view of the fact that there are no free parameters 
(in particular, the overall normalisation is not free) in ei- 
ther curve or data points, and that the assumptions made 
in deriving them are very different. 

Nonetheless, there are systematic differences between 
experiment and theory. The experimental D2 is higher for 
strongly bound particles than its theoretical counterpart by 
a factor up to two. This can be at least partly understood 
from the realisation that the theoretical expression is sin- 
gular at energies close to — ^'(0) and hence fails to describe 
the diffusion process there. The reason for the singularity is 
that in the theoretical derivation there can be no particles 
with E < -*(0) and hence D2{E -*(0)) 0. Yet in 
a finite realisation of a King model, as used in the simula- 
tions, there can be such tightly bound particles and so there 
is no reason why D2 should become zero there. A similar 
argument holds for Di . 

The experimental D's are larger than the theoretical 
ones for particles close to the tidal boundary as well. There 
can be several reasons for that: (1) the initial model is not 
sufficiently close to equilibrium due to the introduction of 
smoothing, and so the effect is purely numerical, (2) the 
effect is real and is due to collective modes or the influence 
of simultaneous encounters, not taken into account in the 
theory, (3) the effect is due to the fact that the effective 
Coulomb logarithm is larger in the outer parts than the value 
used in the theoretical expressions, pmax = Tc- 

Fig. 3 shows the ratios of D2's for various N to D2 
for N = 32k, all for Wo = 9, after taking into account 
the scaling with N suggested by Eq. (7). The shape of the 
experimental D2 curves is very similar for different TV's over 
most of the energy range and in addition they follow the 
expected scaling oc TV/ In A very well, making these ratios 
Si 1. Consequently, all D2's show an increase with respect 
to the theoretical one for loosely bound particles by a factor 
up to 1.5, which suggests that this effect already noted for 
the N = 32k case is real and theory underestimates the 
diffusion rate in the outer parts of this concentrated (Wo = 
9) model. Such an underestimate implies a corresponding 
underestimate of the evaporation rate for this model. 

A comparison between experimental Di's and their the- 
oretical counterparts for various N is made in Fig. 4 from 
which it is clear that these experimental values follow very 
well the theoretical prediction. In particular, there is no ob- 
vious sign of Di being larger in the outer parts, as there was 
in the N = 32k case. Di for core particles in the N = 512 
case fails to track the theoretical curve, unlike the other N 
models, which might be due to the onset of evolution in this 
fewer N system. 

Finally, Fig. 5 compares D's for a different central po- 
tential. Wo = 3 and N = 2k. The correspondence for D2 is 
not as good as in the Wo = 9 case, with the experimental 
D2 typically 50% too low in the outer parts. The correspon- 
dence for Di is better but the statistics are poorer than in 
the Wo = 9 case. 

The diffusion coefficient D2 is a measure of the energy 
relaxation time Te{E) (e.g., Spitzer 1987, Eq. (2-61)), 



Te{E) = {v'{E)}1/D2{E), 



(20) 



where {v^(E))v = 3q/p is the phase-space averaged veloc- 
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(E-Eo) / ^^'(0) 

Figure 1. Relative number of particles at a given energy N[E)/N obtained from Eq. (18) (symbols) compared against M[E)/M, the 
fractional mass at energy E, for a King model (drawn lines) for various and central concentrations Wq , indicated in the figures. Error 
bars assuming Poissonian errors only are indicated in the Wq = 3 figure. Such errors are smaller than the symbols in the other figures. 



ity dispersion of stars with energy E. In view of the good 
correspondence between theoretical and experimental 1)2 's 
we find good agreement between theoretical and experimen- 
tal relaxation rates, although theory slightly underestimates 
Te both in the central parts and in the halo for the con- 
centrated model (factor Si 1.5 — 2). The correspondence is 
slightly worse in the Wo = 3 case. 



6 DISCUSSION 

In the previous section it was shown that there is in general 
excellent agreement between the theoretical and experimen- 
tal diffusion coefficients. In particular, the shape and nor- 
malisation of the Dt(E) curves from both approaches agree 
extremely well (see Fig. 2), using the standard value for the 
Coulomb logarithm. In addition, the experimental D's fol- 
low the theoretically predicted scaling with particle number 
N very well (see Fig. 3). 

However, in spite of this good agreement over most of 
the energy range, the experimental measurement of D2 is 
slightly higher (factor ~ 1.5) in the outer parts of the con- 
centrated model, although no such deviation is apparent for 
Di in those same models (Fig. 4). The measurements of 
D2 in the outer halo for the same King model with differ- 
ent numbers of particles are consistent amongst themselves, 
hence, there seems to be a difference between 'experimen- 
tal' values (based on Eq. (19)) and the 'theoretical' values 
(from Eqs. (6) & (7)). What causes this small discrepancy? 
In the following sections we will elaborate on possible rea- 
sons. Our conclusion from the discussion below is that (1) 
since the definition of states does not involve any dimen- 
sional quantities, it is unlikely to cause a difference between 
core and halo particles and (2) comparison of the models for 
different N and different Wo suggests that the (small) devi- 
ation is not likely to be due to purely numerical effects like 



e.g. mean field relaxation. In addition, the theory involves 
an ad hoc parameter (the Coulomb logarithm) which leads 
us to suggest that the present form of the theory should 
not be regarded as the exact value to compare against, or, 
in other words, the small discrepancy may be an indication 
that both experimental and theoretical estimates of D2 are 
approximations to the real diffusion coefficient. Finally, note 
that the usage of the the Coulomb logarithm does introduce 
length scales, namely po and pmax (section 2.2), so the the- 
ory, in contrast to the experiment, does treat core and halo 
particles on a different footing. 

6.1 Definition of states 

When describing our procedure for measuring the diffusion 
coefficients (section 4.1) it was mentioned that one should 
consider the identification of states with the period between 
two local energy maxima of the single particle as an Ansatz. 
How good is this Ansatzl We will first recall some standard 
results from diffusion theory. 

Suppose a particle undergoes a diffusion process (Brow- 
nian motion) because it undergoes uncorrelated changes in 
position dx over time scales dt. Sampling the position x 
and the associated changes in position Ax of this particle at 
time intervals At, one recovers the standard result that the 
average distance L the particle will wander from its start- 
ing position in a time T grows as = T{Ax^ / At). This 
assumes that position changes of the particle are random 
between successive measurements of its position, i.e., it as- 
sumes that dt « At. However, when the latter inequality is 
satisfied, the value of the diffusion coefficient D = (Ax^ / At) 
is independent of the sampling interval At . 

Next, let us complicate matters and assume that the 
properties of the position change, dx, depend on position x. 
To be more specific, assume these properties change signif- 
icantly when X changes by AL. Consequently, the diffusion 
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N = 32k, Wo=9 Wo=9 




Figure 2. Left scales: comparison between theoretical (curves) 
and experimental (symbols) energy-diffusion rates for the N=32k, 
Wq = 9 case. The diffusion coefficients are scaled so as to 
make them independent of A^. Right scales: idem for the ratio 
of experimental over theoretical D's for loosely bound particles 
[[E — Eo)/^[0) > —0.5) indicated as 'exp. /theory' . Top graph 
refers to D2, lower graph to Di. Error bars are 1(T Poissonian 
errors only. 



coefficient itself will depend on x, D = D{x). It will be clear 
that to measure D{x), one has to sample the position of 
the particle sufficiently often, such that V Ax^ << AL, or, 
At(x) << Td{x} = AL^ / D(x). This expresses the fact that 
the sampling interval At{x) should be much smaller than the 
actual characteristic diffusion time Td{x). Summarising, the 
sampling interval needs to be such that dt{x) << At{x) << 
Td{x). If this inequality is satisfied, the measured value of 
D{x) will not dependent on the actual value At{x) used. 

This example can be applied to energy diffusion in a 
stellar system. However, there is a caveat: which time in- 
terval should be identified with dt, the duration over which 
energy changes are correlated? If most of the relaxation is 
due to near encounters one could presumably take dt ~ X/v, 
where A is the (local) interparticle distance and v the (lo- 
cal) velocity dispersion, because the positions of near par- 
ticles will be uncorrelated when studied intervals of time 
~ dt apart. However, if relaxation is due to collective pro- 



Figure 3. Ratios of diffusion coefficient N D2 / I'a h for given N 
over that quantity for N = 32k and Wq = 9, showing the expected 
scaling cx JV/lnA. Top to bottom: N = 8k, N = 2k, N = 512. 

cesses (e.g. due to interaction of particles with a collective 
oscillation of the system as occurs during violent relaxation), 
energy correlations and hence dt could be much longer, even 
comparable to the local dynamical time ~ l/i/p. 

Not much is known about the time correlation of single 
particle energies in TV-body systems. In making our defini- 
tion of particle states as the time between two successive 
maxima in E{t), we have used the properties of E{t) it- 
self to define when energy changes appear to cease to be 
correlated. / have not proven in this paper that this is in 
fact the case. This assumption is solely made plausible by 
suggesting that after E has reached a local maximum, it 
has 'forgotten' what caused its previous minimum. For ex- 
ample, if a long time scale energy change (an oscillation 
of the system?) dominates the energy change of the parti- 
cle, the states definition will naturally select this long time 
scale for dt. If, however, most of the energy change is due 
to the close passage of another particle, the state will be 
identified as the time over which this other particle has a 
major influence on the particle under consideration. In ad- 
dition to this admittedly hand waving justification of the 
definition, there are some independent points in favour of 
the choice. Firstly, since there are no dimensional quanti- 
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Figure 4. Comparison between experimental and theoretical 
Di 's as in lower panel of Fig. 2 but for different A^'s indicated 
in each panel. 



Figure 5. Same as Fig. 2 but for Wq = 3 and showing ratios for 
loosely bound particles with {E - Eo)/'^(0) > -0.4. 



ties involved, the definition treats all particles (e.g. in core 
or halo) in exactly the same manner. Consequently, even if 
the inequality dt << At were nof satisfied, it would appear 
mysterious that the resulting experimental 1)2 's fit the the- 
oretical ones in the core, but not in the halo. In addition, 
choosing too small values of At can only lead to an overes- 
timate of D2, since energy changes are added in quadrature 
to compute this quantity. Consequently, such an argument 
does not allow one to explain why the measured values of 
D2 appear too small compared with the theoretical values 
for the Wo = 3 model. Secondly, from the measured values 
of D2 it is clear that the second equality. At << Te is al- 
ways satisfied: in fact, if it were not, the process could not 
be considered a diffusion process in the first place! (Note 
in passing that it is impossible to satisfy both inequalities 
with a single (energy independent) time step: indeed, the 
core relaxation time Ti;(core) is actually shorter than the 
characteristic time scale X/v in the halo, making it impos- 
sible for a single At to satisfy dt{E) « At « Te(E) 
throughout the system.) In addition, note that our defini- 
tion of states allows us to measure Di{E) as well. To our 
knowledge it is the first time that this quantity has been 
measured directly, and its measured value agrees excellently 



with the theoretical prediction for all N and Wo , both in 
core and halo. 

Finally, as we remarked before, the theoretical expres- 
sions are singular, both in the very core of the system 
(E *(0)) and in the outer halo (E Eo). This is be- 
cause the theory does not take into account that in an actual 
TV-body realisation of the King model, there can be parti- 
cles more strongly bound than ^'(0) or more loosely bound 
than Eo , unlike in the analytic model. This singular nature 
of the theoretical D's causes the discrepancy most clearly 
seen in e.g. Fig. 3 for E ^(0)- Clearly, at least in this 
case, it is the theoretical expression which is to blame for 
the discrepancy! 

We conclude that, although the assumption that energy 
changes are uncorrelated between successive states (where 
states are defined in section 4.1) remains an Ansatz until 
further work is done, the fact that the definition treats halo 
and core particles in exactly the same way suggests that 
the small discrepancy between experiment and theory in the 
halo of the models is unlikely to be due to a failure of the 
present method. 
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6.2 Simulations 

The experimental curves compare the diffusion coefficients 
of the numerical King model with those of an analytic King 
model of specified concentration. Is this comparison fair, i.e., 
is the numerical King model a good representation of the an- 
alytic model? In turn, we will discuss the importance of dis- 
creteness {l/\/W noise), possible out-of-equilibrium effects 
(breathing modes) and the change in structure as a conse- 
quence of the relaxation process itself (leaving the conse- 
quences of the rate of change of structure due to relaxation 
effects to the next section). 

When generating a King model with a finite number 
N of particles one naturally suffers from I/t/TV noise and 
one would expect the model to be a better representation 
of the analytic system for larger values of N. However, as 
Fig. 3 testifies, the experimental values of D2 for different 
N are consistent amongst themselves, over the range 64 in 
N shown. In addition, such l/\/W noise would be unable to 
explain why experimental 1)2 's fit in the core (which has few 
particles) and not in the halo. 

A numerical model generated using a random number 
generator is never completely in equilibrium. In the present 
case, this is aggravated by the fact that the analytical model 
does not take into account the smoothing employed in the 
numerical model. Consequently, it might be expected that 
the numerical King models do not start in equilibrium. How- 
ever, as in the previous case, the amount by which this influ- 
ences the results should depend on N, since the numerical 
model actually converges to the analytic one for N ^ 00. 
Consequently, the amount that such breathing modes con- 
tribute to the measured value of D2 should be different for 
different N, yet no such effect is seen in Fig. 3: it appears 
that the non-equilibrium initial state from which these mod- 
els are started does not overly influence the measurement of 

Finally, we have compared numerically evolved models 
with a specific unevolved King model. But how much has 
the relaxation process itself changed the numerical models? 
The half mass relaxation time TE{h) is known to be of or- 
der N/ In A times the dynamical time td (e.g. Spitzer (1987), 
BT87), hence, f» 10* for the N = 32k model, using the 
measured value In A Si 3 for the Coulomb logarithm. Con- 
sequently, we have evolved that model over only 0.1% of 
its relaxation time: clearly, relaxation itself is completely 
unimportant for the N = 32k case and is not able to ex- 
plain the small discrepancy in the halo where the relaxation 
time is likely to be even larger than TE{h). The unimpor- 
tance of relaxation itself is born out by the behaviour of 
the Lagrangian radii for this model, which remain virtually 
identical to the analytic ones over the whole time-span the 
simulation lasted. In addition, since the relaxation time is 
a strong function of N , the simulated time for the different 
models, in units of their relaxation time, is very different, 
e.g., the N = 512 model has been integrated over ~20% of 
its relaxation time (so 200 times longer than the N = 32k 
model), yet the measured value of D2 agrees well with the 
measurement from the N = 32768 case in the halo. 



6.3 Mean field relaxation 

The rate at which the single particle energy E{t) changes 
is a sum of the rate due to diffusion, as quantified by the 
theoretical expression for Di , and due to mean field relax- 
ation, which changes the potential as a function of radius 
and hence E{t) as well. It is the sum of these two that is 
measured by the experimental Di. Specifically, (dE/dt)T = 
{dE/dt}MF + {dE/dt}D (T for total, MF for mean field, D for 
diffusion) and Eq. (2) is a measure of (dE/dt)D but Eq. (19) 
measures (dE/dt)T- An estimate for [dE I dt) mf can be ob- 
tained using the Fokker-Planck simulations of Cohn (1979, 
1980), who gives an estimate of the logarithmic rate of 
change of the central density, /9(0), in units of the central 
relaxation time tr(0): 



Using Figs. 2 and 3 in Cohn (1979) to convert the central 
relaxation time to the half-mass relaxation time tru = F x 
tr{0), with F w 100 for a Wo = 9 King model, we estimate 
for the rate of change in the centre: 

l/'^^ro^N ^''M ^ | <i^(0) tru |_| '^ln(^(0)) , I 

« \^M£mtr,\=Fi^0.5 (22) 

using the value { f» 5 x 10"^ (Fig. 6 in Cohn 1980); trh is 
the initial half-mass relaxation time, which corresponds to 
a Plummer model and so is likely to be larger than for the 
concentrated King model. (Cohn's simulations start from a 
Plummer model but he shows that the later stages of evo- 
lution can be fit reasonable well by members of the family 
of King models.) In Fig. 6 we plot (dE/dt)T = E>i{E) ob- 
tained from Fig. 2 and converting the dynamical time to 
half-mass relaxation time using tru = Ntd/8lnA (BT87, 
Eq. 4-9). Comparing the numerical value in Eq. (22) with 
Fig. 6 it is clear that the central value \(dE /dt{0)) mf\ is 
at least an order of magnitude smaller than the central 
value of \(dE /dt)T\- In addition, the large difference be- 
tween (dE /dt{0)) MF and (dE/dt)T suggest that (dE/dt)T ~ 
[dE I dt) B , which explains the good correspondence in Fig. 4. 
We conclude that the contribution to (dE/dt)T of the mean 
field relaxation rate is likely to be a small fraction (< 1% 
say) of the contribution due to diffusion. 

We conclude from these remarks that, although the 
use of states to measure the diffusion coefficients is based 
on an unproven assumption, there appear to be no serious 
shortcomings in either method or numerical simulations that 
would easily explain why the experimentally measured dif- 
fusion coefficient should be off in the halo, even though they 
fit well in the core. We re-iterate however, that we consider 
the major conclusion to be drawn from Fig. 2 is that there 
is excellent agreement between experimental and theoretical 
diffusion coefficients. More work is needed to judge whether 
the factor Si 1.5 disagreement in the outer halo is due to 
hidden inaccuracies in the present measurement of D2 or is 
due to a slight underestimate by the standard theory. In this 
connection we suggest that the fact that one has to introduce 
ad hoc a cut-off in an integral appearing in the theoretical 
derivation (the Coulomb logarithm) suggests that improve- 
ments in the theory might be possible. For example, if one 
introduces a cut-off for distant encounters, it would appear 
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Figure 6. Experimental value for {dE/dt)x = Di{E) from 
Eq. (18) and (19) (symbols) against the theoretical expression 
for {dE / dt) jj obtained from Eq. (2) (drawn lines) for the concen- 
trated Wq = 9 King model with = 32k particles. As before, 
errorbars denote Poissonian errors only. The inset zooms in on 
the more loosely bound part of the system. 



evident that such a cut-off scale should naturally depend on 
the location of the particle. We also noted in the Introduc- 
tion that the very concept of a single Coulomb logarithm 
becomes suspect when dealing with a very inhomogeneous 
system, as is the case for these concentrated King models. 



7 SUMMARY 

Energy diffusion coefficients based on phase-space averages 
of energy gains per unit time for independent two-body en- 
counters are good approximations to the experimental val- 
ues obtained from TV-body simulations. These experimental 
values are based on the notion of particle states, defined as 
follows: particle i is considered to be in a given state dur- 
ing the time interval between two consecutive local maxima 
Et{tb) and Et{te) of its energy vj l2 — 4>t. The state is then 
characterised by its energy E,{tb), duration te — tb and tran- 
sition energy E,{te) — E,{tb). It is assumed that different 
states are independent for all tb and all particles i. Statis- 
tics on particle states are then used to compute experimental 
diffusion coefficients. Although the agreement between these 
experimental coefficients and their theoretical counterparts 
is good, theoretical diffusion coefficients are systematically 
too small by a factor Si 1.5 — 2 both in the core and in the 
outer halo for the concentrated model, whereas in the less 
concentrated model theory overestimates the diffusion rate 
by a similar factor. 

Simulations that study the dependence on angular mo- 
mentum and the effect of a mass-spectrum will be done in 
the near future. It is important to redo this type of calcu- 
lation using regularization instead of numerical smoothing 
to determine whether the results presented here are quali- 
tatively applicable to globular clusters. 
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